clear 
est clear

********************************************************************************
*	I - Explore attendance and such
********************************************************************************

*** Prep compliance data  *********************

use "${data}\CashTherapy_log.dta", clear
*** one of the CRs was a mentor??? 
drop if CR_ismentor == 1 
* clean up attendance
replace attended_atleast1=. if attended_atleast1<0
label var attended_atleast1 "Therapy, attended at least one"
* plot attendance
tab tot_attendance attended_atleast1 
* hist tot_attendance , d
* hist tot_attendance if attended_atleast1 == 1 , d
* clean up cash 
label var CT_status "Received cash"
rename CT_status cash 
fre cash
* save 
keep cr_id attended_atleast1 tot_attendance cash Week14attendance
tempfile compliance 
save `compliance', replace

*** Prep data with compliance *********************

use "${data}/Analysis_wide.dta", clear 
distinct block club
*** merge compliance, 
*	just drop the non-mergers, won't make a difference as their outcomes are missing anyways...
merge m:1 cr_id using `compliance', nogen assert(match master) keep(match) 
* check
tab attended_atleast1 , m
tab cash, m
*** clean up the label
label list treat_ther
label define treat_ther 1"IPT-G", modify


*** Check out some stuff *********************

*** plot attendance
* cdf
*distplot tot_attendance, over(treat) 
distplot tot_attendance, ///
	graphregion(color(white)) legend(rows(1) pos(6)) lc(green)  /// 
	title("No. of therapy sessions attended.") 
graph export "${figures}/attendance_cdf.png", replace 
* pdf
*hist tot_attendance if treat!=0, d by(treat) 
hist tot_attendance if treat!=0, d ///
	graphregion(color(white)) legend(rows(1) pos(6)) color(midblue) lcolor(grey) lwidth(vthin)  ///
	title("No. of therapy sessions attended.") 
graph export "${figures}/attendance_pdf.png", replace 

twoway ///
    (histogram tot_attendance if treat == 1, discrete start(0) color(red%50) lcolor(none) percent) ///
    (histogram tot_attendance if treat == 2, discrete start(0) color(blue%50) lcolor(none) percent), ///
    legend(order(1 "IPT-G" 2 "IPT-G+") position(6) region(lwidth(none) margin(none))) graphregion(color(white)) ylabel(, angle(0)) xlabel(, angle(0)) ///
    xtitle("No. of therapy sessions attended.") ytitle("Percentage") 
graph export "${figures}/attendance_pdf_treat.png", replace 


 
*** now, check different levels of attendance 
* dummies
gen attended_atleasthalf 	= tot_attendance >= 7	if !missing(tot_attendance)
gen attended_atleast10 		= tot_attendance >= 10	if !missing(tot_attendance)
gen attended_atleast12 		= tot_attendance >= 12	if !missing(tot_attendance)
gen attended_all	 		= tot_attendance == 14	if !missing(tot_attendance)
mean attended* , over(treat)
* make it a categorical 
egen attended_group = rowtotal(attended_atleast1 attended_atleast10 attended_all) 
replace attended_group = -1 if treat == 0 
label define attended_group -1 "Control group" 0 "Attended 0 sessions" 1 "Attended 1-9 sessions" 2 "Attended 10-13 sessions" 3 "Attended all 14 sessions" 
label val attended_group attended_group
tab attended_group

*** mean by attendance groups 
tabstat phq8_min1 phq8_score1 ghq12_min1 ghq12_score1 , by(attended_group)  stat(mean) col(variables) nototal format(%9.3f)
tabstat phq8_min2 phq8_score2 ghq12_min2 ghq12_score2 , by(attended_group)  stat(mean) col(variables) nototal format(%9.3f)
tabstat phq8_min3 phq8_score3 ghq12_min3 ghq12_score3 , by(attended_group)  stat(mean) col(variables) nototal format(%9.3f)

tabout  attended_group using "${tables}/attendance_mean_outcomes1.tex" , sum c(mean phq8_min1  mean phq8_score1 mean ghq12_min1 mean ghq12_score1 N phq8_min1) replace style(tex) bt format(3) ptotal(none) h3(nil) h2(nil)
*tabout  attended_group using "${tables}/attendance_mean_outcomes2.tex" , sum c(mean phq8_min2  mean phq8_score2 mean ghq12_min2 mean ghq12_score2 N phq8_min1) replace style(tex) bt format(3) ptotal(none)
*tabout  attended_group using "${tables}/attendance_mean_outcomes3.tex" , sum c(mean phq8_min3  mean phq8_score3 mean ghq12_min3 mean ghq12_score3 N phq8_min1) replace style(tex) bt format(3) ptotal(none)





 




********************************************************************************
*	II - Impacts of Therapy vs. Control - LATE Estimates
********************************************************************************

*** Prep data  *********************

use "${data}\CashTherapy_log.dta", clear
* clean up attendance
replace attended_atleast1=. if attended_atleast1<0
* label outcomes 
label var attended_atleast1 "Therapy, attended at least one"
rename attended_atleast1 therapy 
fre therapy
label var CT_status "Received cash"
rename CT_status cash 
fre cash
* save 
keep cr_id therapy cash
tempfile compliance 
save `compliance', replace
 
use "${data}/Analysis_wide.dta", clear 
distinct block club
*** is it easier in long mode? I think it might be...
reshape long phq8_score ghq12_score rosb_score res_score loc_score ghq12_min phq8_min w, i(block cr_id club treat ) j(round)
*** drop +cash folks in longer follow ups 
drop if treat==2 & round>=2
tab round 
*** merge compliance, 
*	just drop the non-mergers, won't make a difference as their outcomes are missing anyways...
merge m:1 cr_id using `compliance', nogen assert(match master) keep(match) 
*still one obs missing the therapy, so drop
tab therapy , m
drop if missing(therapy)
tab cash, m
*** clean up the label
label list treat_ther
label define treat_ther 1"IPT-G", modify
label define ment 0 "No" 1 "Therapy", modify


local treatment  	therapy 
local instrument 	treat_ther


*** the main four in all 3 waves 
qui forvalues r = 1/3 {
	foreach y in phq8_min phq8_score ghq12_min ghq12_score {
		*center the blocks within round and non-missing sample
		capture drop blocks* 
		qui forvalues b=1/16 {
			gen blocks`b' = block == `b' 				if !missing(`v') & round == `r' 
			su blocks`b'				 				if !missing(`v') & round == `r'
			replace blocks`b' = blocks`b' - r(mean)		if !missing(`v') & round == `r'
			su blocks`b'				 				if !missing(`v') & round == `r'
		}
		local ctrls blocks1-blocks15
		*** regress
		eststo `y'_`r':	ivreg2 `y' `ctrls' (b0.`treatment' b0.`treatment'#c.(`ctrls') = b0.`instrument' b0.`instrument'#c.(`ctrls')) [pw=w] if round==`r', cluster(club) partial(`ctrls') 
		*eststo `y'_`r':	ivreg2 `y' `ctrls' (b0.`treatment' = b0.`instrument') [pw=w] if round==`r', cluster(club) // partial(`ctrls') 
			*** add scalars
			* control mean
			qui su `y' [aw=w] if `instrument' == 0 & round==`r' 
			estadd scalar c_mean = r(mean) 	
			estadd scalar c_sd = r(sd)
			* p-value 
			test 1.`treatment' == 0 
			estadd scalar pval = r(p) 
			* first-stage f-stat 
			estadd scalar ffstat = e(cdf) 
	}
	drop blocks*
}
*** the secondary 4 in waves 1 and 3
qui forvalues r = 1(2)3 {
	foreach y in rosb_score res_score loc_score {
		*center the blocks within round and non-missing sample
		capture drop blocks* 
		qui forvalues b=1/16 {
			gen blocks`b' = block == `b' 				if !missing(`v') & round == `r' 
			su blocks`b'				 				if !missing(`v') & round == `r'
			replace blocks`b' = blocks`b' - r(mean)		if !missing(`v') & round == `r'
			su blocks`b'				 				if !missing(`v') & round == `r'
		}
		local ctrls blocks1-blocks15
		*** regress
		eststo `y'_`r':	ivreg2 `y' `ctrls' (b0.`treatment' b0.`treatment'#c.(`ctrls') = b0.`instrument' b0.`instrument'#c.(`ctrls')) [pw=w] if round==`r', cluster(club) partial(`ctrls') 
		*eststo `y'_`r':	ivreg2 `y' `ctrls' (b0.`treatment' = b0.`instrument') [pw=w] if round==`r', cluster(club) // partial(`ctrls') 
			*** add scalars
			* control mean
			qui su `y' [aw=w] if `instrument' == 0 & round==`r' 
			estadd scalar c_mean = r(mean) 	
			estadd scalar c_sd = r(sd)
			* p-value 
			test 1.`treatment' == 0 
			estadd scalar pval = r(p) 
			* first-stage f-stat 
			estadd scalar ffstat = e(cdf) 
	}
	drop blocks*
}

* esttab *1 , b se keep(1.therapy) stat(c_mean N ffstat pval_ar) star(* 0.10 ** 0.05 *** 0.01)
* esttab *2 , b se keep(1.therapy) stat(c_mean N ffstat pval_ar) star(* 0.10 ** 0.05 *** 0.01)
* esttab *3 , b se keep(1.therapy) stat(c_mean N ffstat pval_ar) star(* 0.10 ** 0.05 *** 0.01)

*** unadjusted table
* RR 
#delimit ;
esttab 	phq8_min_1 ghq12_min_1 
		phq8_score_1 ghq12_score_1 	 
		rosb_score_1 res_score_1 loc_score_1
	using "${tables}/results-mh-late1.tex", replace  booktabs
	label nostar nonote nomtitle nonumber frag nolines nogaps
	b(3) se(3) keep(1.therapy) 
	stats(c_mean c_sd N pval ffstat, fmt(3 3 0 3 3)  
	label("Control mean" "Control SD" "Observations" "H0: IPT-G=0 p-value" "First-stage F-stat.")) 
	; 
#delimit cr 
* Midline 
#delimit ;
esttab 	phq8_min_2   ghq12_min_2
		phq8_score_2 ghq12_score_2 	
	using "${tables}/results-mh-late2.tex", replace  booktabs
	label nostar nonote nomtitle nonumber frag nolines nogaps
	b(3) se(3) keep(1.therapy) 
	stats(c_mean c_sd N pval ffstat, fmt(3 3 0 3 3)  
	label("Control mean" "Control SD" "Observations" "H0: IPT-G=0 p-value" "First-stage F-stat.")) 
	; 
#delimit cr 
* Endline 
#delimit ;
esttab 	phq8_min_3  ghq12_min_3  	
		phq8_score_3 ghq12_score_3 	 
		rosb_score_3 res_score_3 loc_score_3
	using "${tables}/results-mh-late3.tex", replace  booktabs
	label nostar nonote nomtitle nonumber frag nolines nogaps
	b(3) se(3) keep(1.therapy) 
	stats(c_mean c_sd N pval ffstat, fmt(3 3 0 3 3)  
	label("Control mean" "Control SD" "Observations" "H0: IPT-G=0 p-value" "First-stage F-stat.")) 
	; 
#delimit cr 
*/
